Orbital magnetoelectric coupling at finite electric field 
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We extend the band theory of linear orbital magnetoelectric coupling to treat crystals under finite 
electric fields. Previous work established that the orbital magnetoelectric response of a generic 
insulator at zero field comprises three contributions that were denoted as local circulation, itinerant 
circulation, and Chern-Simons. We find that the expression for each of them is modified by the 
presence of a dc electric field. Remarkably, the sum of the three correction terms vanishes, so that 
the total coupling is still given by the same formula as at zero field. This conclusion is confirmed 
by numerical tests on a tight-binding model, for which we calculate the field-induced change in the 
linear magnetoelectric coefficient. 

PACS numbers: 75.85.+t,03.65.Vf,71.20.Ps 



c 

o 

(N 
> 

NO 
OO 
NO 
^1" 



X 



Magnetoelectrics are magnetic insulators whose dielec- 
tric polarization P changes linearly under a small applied 
magnetic field B and, conversely, whose magnetization 
M changes linearly with a small applied electric field 
£. 1,2 This linear magnetoelectric (ME) coupling is de- 
scribed by the response tensor 3 

which is odd under both spatial inversion [V) and time- 
reversal (T) symmetries. Thus ME materials must be 
acentric and display magnetic order. 

In crystals where only one of the two symmetries, V or 
T, is present, it may still be possible to induce a linear 
ME effect by applying an external field which breaks that 
symmetry. So, for example, a centrosymmetric insulating 
antiferromagnet placed in a (strong) electric field loses its 
inversion center. Likewise, a nonmagnetic ferroelectric 
crystal loses time-reversal symmetry when subject to a 
magnetic field. In both cases the symmetry is sufficiently 
lowered that the tensor a becomes nonzero. 

It is useful to view these field-induced effects as higher- 
order ME responses of the unbiased crystal. 4 Two quad- 
ratic ME effects can be defined in this way. Going to 
next order in magnetic field yields the tensor 



a _ fay 

Pijh ~ 8B k 



dBidBk : 



(2) 



which is odd under V and even under T. Going instead 
to next order in the electric field gives 



lijk 



d£ k 



d 2 Mj 
d£jd£ k 
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which is even under V and odd under T. Reference 4 
lists the form of these tensors for all the crystal classes. 



While most investigations of ME couplings in solids have 
focused on the linear response a for a reference state 
of the crystal at zero electric and magnetic fields, the 
quadratic responses /? and 7 have also been measured 
in materials where a vanishes by symmetry. In particu- 
lar the electric-field- induced effect, which constitutes the 
primary focus of this work, was first measured by O'Dell 
in yttrium iron garnet. 5 

The ME response can be divided into four contribu- 
tions, depending on whether the response is frozen-ion 
(purely electronic) or lattice-mediated, and whether it is 
spin or orbital in character. We will refer to the frozen- 
ion part of the orbital response as the orbital magneto- 
electric polarizability (OMP). 6 ' 7 While the OMP is typ- 
ically a small contribution to the ME response in con- 
ventional magnetoelectrics, it was recently realized that, 
under certain conditions of surface preparation, Z2-odd 
topological insulators 8 should display a large, quantized 
OMP response. 6 ' 9 This is a remarkable prediction, espe- 
cially considering that in this class of materials T symme- 
try is preserved in the bulk (it must, however, be broken 
on the surface). This topological magnetoelectric effect 
has triggered a great deal of interest in orbital magneto- 
electric couplings in solids. 

The microscopic theory needed to calculate the OMP 
at zero electric and magnetic fields from first principles 
was worked out in Refs. 7 and 10. In addition to the so- 
called Chern-Simons term responsible for the topological 
ME effect, 6,9,11 it was found that two more (Kubo) terms 
contribute to the OMP in conventional magnetoelectrics 
in which T and V symmetries are broken spontaneously 
in the bulk. 

In this work we generalize the band theory of OMP of 
periodic insulators 7,10 to finite electric fields. That is, we 
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evaluate the coefficient a at nonzero £, 



d£,. 



(4) 
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(Henceforth, the condition B = will be implied 
throughout. It is also understood that from now on a de- 
notes the OMP part of the entire ME response.) A prin- 
cipal result of our work is the conclusion that the zero- 
field expression for the total OMP remains valid at finite 
electric field, while the above-mentioned Chern-Simons 
and Kubo terms separately acquire field-induced contri- 
butions. We confirm our formal results by numerical tests 
on a tight-binding model. 

Our derivation of a formula for a(£) proceeds along 
the lines of Ref. 7. We start from the expression given 
therein for the orbital magnetization of a generic band 
insulator under a finite electrical bias. It comprises three 
terms, 



M 3 {£) = Mf G (£) + M iu (£) 



IC / 



(£), (5) 



where 



ikT LC = 
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d 3 klm(d p u nk \H^.\d q u nk ), (6) 



and 



d 3 klm j(w„k|#kl u mk)(dpW m k|dgU„k) j 

(7) 



M CS = e J]_ £ / d 3 k tj . 
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The common prefactor in these formulas is rj = 
— e/S(27r) 3 (e > is the magnitude of the electron 
charge), and a sum is implied over repeated Cartesian 
(pqr) and valence-band (mn) indices. The cell-periodic 
part of the field-polarized Bloch state 12 is denoted by 
\u n u), dj is the partial derivative with respect to the 
jth component of the wavevector k, and the tilde in- 
dicates a covariant derivative dj = Qkdj, where Qk = 
1 — |«nk)(«nk| (sum implied over n). The Hamiltonian 
is defined as 



-ikr^Ogik-r 



(9) 



where H° is the zero-field part of the crystal Hamiltonian. 
In Eq. (8) the symbol A p denotes the Berry connection 
matrix 



on the electric field is only implicit, M cs displays an ex- 
plicit linear dependence on £. It is therefore expedient 
to introduce a new quantity M-p s via the relation 



Mf 3 {£) =£jM^ s (£), 



(11) 



where the subscript '1' serves as a reminder that M GS 
enters the expression for M multiplied by £ to the first 
power. 

All three magnetization terms, M LC , M IC , and M cs , 
are invariant under gauge transformations within the 
valence-band manifold, although in the case of M cs this 
invariance is only modulo a quantum of indeterminacy. 9 
In the limit that £ goes to zero, M cs vanishes and Eq. (5) 
reduces to the expression for the spontaneous orbital 
magnetization. 13 

As already mentioned, all terms in Eq. (5) can con- 
tribute to the linear ME coupling, Eq. (4), so that 



a ij (£) = af 1 c (£) + a%(£) + a??(£). 



(12) 



The derivation of the expressions for these objects is 
straightforward though somewhat lengthy. It essentially 
repeats the steps in Appendix B of Ref. 7, where the 
derivation was carried out for the LC and IC ("Kubo") 
terms under the assumption that £ = (the CS term is 
trivial at £ = 0). At £ ^ one may show that each of 
the terms in Eq. (12) consists of a "zero-field" part plus 
a "field-correction" part having an explicit linear depen- 
dence on £, 



un{£) = a .ij{£) + £ja hi (£). 



(13) 



The field-correction terms for the LC and IC contribu- 
tions can be traced back to Eqs. (B.7) and (B.8) in Ref. 7, 
which at £ ^ acquire extra terms. As for the Chern- 
Simons contribution, differentiating Eq. (11) with respect 
to £j yields a%% = SijMf 8 and a°f = dM GS /d£ t . 
Thus, we arrive at the results 



*0%(£) = VCjpql™ j rf 3 fc((^nk|(9,^)|Ali„k) 

- 2 1 (^P u nk|(A#k)l^ U ™k)), 



(14) 



a 0?ij(£) = ^pg Im J d3k ((3pU„k|AWmk}(Wmk|(dg#k)Kk 
~\ 0pU n k\d q U mk ) (w mk I ( A#k) \Unk)) ) 



(15) 
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and the trace is over the valence bands. 

Equations (6) and (7) describe respectively the 
local and itinerant circulation contributions to the 
magnetiztion, 7 while Eq. (8) is the Chern-Simons term. 
At variance with the other two terms, whose dependence 
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and 
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In the above expressions, Di is the partial derivative with 
respect to the ith component of the electric field. The 
terms containing DiH^. in Eqs. (14) and (15) are screen- 
ing corrections which are present in self-consistent calcu- 
lations. 

Equations (14)-(16) for the zero-field terms are essen- 
tially rewritten from Ref. 7. It should be emphasized, 
however, that in the present context these expressions 
depend on the electric field implicitly via the wave func- 
tions. The explicit field dependence is given by the field- 
correction terms, Eqs. (17) and (18). Remarkably, these 
terms are not independent and add up to zero when in- 
serted into Eq. (12). We conclude, therefore, that the 
expression for the total OMP derived in Refs. 7 and 10 
assuming £ = remains valid for £ ^ 0. This consti- 
tutes one of our principal results. The explicit expres- 
sion given in Eq. (17) for the field-correction terms is 
the other main result of this work. It is useful if one is 
interested in the field dependence of the separate gauge- 
invariant contributions to the OMP. Because it contains 
three k derivatives and one field derivative, this quantity 
is even under V and odd under T, just like the coeffi- 
cient 7 defined by Eq. (3). This is reasonable since, as 



one can see from Eq. (13), a 1 



LC/IC/CS 



gives a contribu- 



tion to ^ LC / IC / CS anc i should therefore have the same 
symmetry properties. 

As a check of our analytic derivation, we have imple- 
mented the formula for a{£) in a tight-binding model, 
and used it to calculate the nonlinear ME coefficient j zzz 
at £ = 0. Since the tensor 7 vanishes in T- invariant sys- 
tems, we need a model where T is spontaneously broken, 
and we chose that of Ref. 7. This is a spinless model 
with eight sites per primitive cell arranged on a 2 x 2 x 2 
cube, where T symmetry is broken by complex nearest- 
neighbor hoppings, and we have used the same on-site 
energies and nearest-neighbor hoppings tabulated in that 
work. (This choice of parameters also breaks V, so that 
the linear ME tensor a is nonzero already at £ = 0, but 
this is not essential for our present purposes.) As in Ref. 7 
the two lowest bands were treated as occupied, and the 
phase (p of one of the complex hoppings was chosen as a 
control parameter for plotting purposes. 

The technical details of the tight-binding implementa- 
tion of Eqs. (6)-(8) and (14)-(17) can be found in Ref. 7. 
The only significant difference with respect to that work 
is that the field derivative |-DiU n k) of the cell-periodic 
Bloch states must be evaluated at finite £. Under these 
circumstances the usual "sum-over-states" formula 13 can- 
not be employed, and one must instead minimize a suit- 
ably defined functional. 14 

We shall calculate the zzz component of 7 from the 
first equality in Eq. (3). Combining with Eq. (12) we 
find 



_, _,LC 1 _,IC 1 _,CS 

7=7 +7 +7 . 



(19) 
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FIG. 1. Decomposition of 7 ZZZ of Eq. (19) into j LC (solid 
lines), 7 IC (dashed lines), and 7° s (dotted lines) calculated 
using Eqs. (20) and (21). Symbols denote the same contribu- 
tions evaluated using Eq. (22). 



The CS term is the simplest to evaluate, as the deriva- 
tive of Eq. (13) with respect to £ z can be taken analyt- 
ically. The zero-field and field-correction terms therein 



both contribute an amount z (0) to 7^(0). Thus 



7^(0) = 2af?(0) = -4a£(0), 



(20) 



where the second equality follows from Eq. (18). The 
quantity on the right-hand side can be evaluated directly 
from Eq. (17). For the LC and IC terms we calculate the 
derivative of the zero-field terms in Eq. (13) using finite 
differences and obtain 



LC/IC/~ 



-LC/IC 

I ZZZ 



(0) 



LC/IC 



+ a\ c z (0). (21) 



2£ z 



In practice we evaluate the first term from Eqs. (14) and 
(15), using small positive and negative fields along z of 
magnitude £ z = 1.0 x 10~ 5 V/m. 

The results of the above calculations were compared 
with a finite-difference determination of the second field 
derivative of M, 



1zzz(0) 



d 2 M z 



d£l 



£=0 



M z (£ z ) - 2M Z (0) + M z (-£ z ) 
£1 



(22) 

using the fc-space expressions from Ref. 7 for the LC, IC 
and CS terms in Eq. (5). The results obtained in this 
manner can be taken reference, since the fc-space 
expression for M(£) has been carefully tested by com- 
paring with real-space calculations on bounded samples 
cut from the bulk crystal. 7 

The agreement between the two sets of calculations can 
be seen in Fig. 1, where the LC, IC, and CS contributions 
to 7z2z are plotted separately as functions of <p. In this 
calculation ^ zz is about an order of magnitude smaller 
than 7^.. From Eqs. (20) and (21) it then follows that 
the field-correction terms contribute little, especially in 
the case of 7^. Further numerical tests focusing on 
those small terms are therefore desirable. 
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FIG. 2. (a) Right-hand side (solid line) and left-hand side 
(symbols) of Eq. (23) . Squares and circles denote the LC and 
IC contributions, respectively, (b) Equation (20) (solid line) 
and Eq. (22) for the CS contribution (crosses), both multiplied 
by a factor of —1/4 for visual check of Eq. (18) by comparison 
to (a). 



In order to isolate the field-correction terms in 7^ 
and 7^22, we subtract the zero- field terms from the total: 



da 



LC/IC 

0,zz 



d£ z 



LC/IC 



(0). (23) 



In Fig. 2 (a) we plot, as a function of ip, the two sides of 
this equation. The field derivatives on the left-hand side 
are evaluated by finite differences, while the right-hand 
side is calculated from Eq. (17). It is clear that the field- 
correction terms in Eq. (13) are nonzero, and the good 
agreement between the three curves demonstrates that 
for both LC and IC they are given by Eq. (17). 

The CS contribution does not need additional tests 
since, as noted above, the contributions to 7^ from the 
zero-field and field-correction terms are identical. How- 
ever, we reproduce in Fig. 2 (b) the CS curve from Fig. 1 
multiplied by a factor —1/4, so that the correctness of 
Eq. (18) can be verified by direct visual inspection. This 
completes the numerical checks of the A:-space formula 
for a(£). 

To summarize, we have extended the recently devel- 
oped band theory of orbital magnetoelectric response to 
treat crystals under a finite electrical bias. The theory 
presented in this work may be especially useful in calcu- 
lations of the second-order magnetoelectric effect defined 
by Eq. (3). While it is possible in principle to calculate 
the second derivative of M by finite differences, the nu- 
merical stability is likely to be improved by taking one of 
the field derivatives analytically, leaving only one deriva- 
tive to be performed numerically We have demonstrated 
that in order to calculate the total OMP at finite electric 
field, one may use the same equations (14)-(16) that were 
previously derived for zero field. This is true even though 
the individual local-circulation, itinerant-circulation, and 
Chern-Simons contributions do separately acquire field- 
correction terms. At present, we are not aware of any 
simple argument that could have anticipated the exact 
cancellation of these terms in the expression for the total 
OMP. 
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